Large effect of polydispersity on defect concentrations in colloidal crystals 
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We compute the equilibrium concentration of stacking faults and point defects in polydisperse hard-sphere 
crystals. We find that, while the concentration of stacking faults remains similar to that of monodisperse hard 
sphere crystals, the concentration of vacancies decreases by about a factor two. Most strikingly, the concentra- 
tion of interstitials in the maximally polydisperse crystal may be some six orders of magnitude larger than in a 
monodisperse crystal. We show that this dramatic increase in interstitial concentration is due to the increased 
probability of finding small particles and that the small-particle tail of the particle size distribution is crucial for 
the interstitial concentration in a colloidal crystal. 



I. INTRODUCTION 

The experimental study of colloidal crystals is of interest 
for at least two reasons. First of all, the possibility to design 
the constituents of such crystals, allows us to gain insight into 
the factors that determine the structure and kinetics of forma- 
tion of crystalline materials. In addition, colloidal crystals are 
of interest because of their potential application as photonic 
materialsl 1|. To a first approximation, one might view col- 
loidal crystals as scale models of atomic crystals. But this 
analogy is flawed for several reasons. First of all, the inter- 
molecular forces between colloidal particles may be qualita- 
tively different from those between atoms. Secondly, the dy- 
namics of colloidal matter is intrinsically different from that of 
atomic materials, due to the presence of a solvent. Finally, un- 
like atomic materials, colloidal systems are never completely 
monodisperse. This polydispersity may have important con- 
sequence for the phase behavior and structural properties of 
the colloidal crystals. In addition, polydispersity can have an 
effect on the equilibrium concentration of (point) defects in 
colloidal crystals. As defects may strongly influence the pho- 
tonic properties of colloidal crystals, a better understanding of 
the effect of polydispersity on defect concentrations, may also 
be of practical relevance for the design of photonic crystals. 

In the present paper, we describe a numerical study of the 
effect of polydispersity on the concentration of stacking faults, 
vacancies and interstitials in hard-sphere colloidal crystals. 



II. SIMULATION METHODS 

A. Semigrand Canonical Ensemble 

To simulate a the equilibrium properties of polydisperse 
hard-sphere cr^tals, we used the semigrand canonical ensem- 
ble method |2, 3]. For a system with continuous size polydis- 
persity, the free-energy functional of the semigrand canonical 
ensemble is given by: 

Y{N,PJ,aQ,{Ali]) = U-TS + PV+Nii (Oo) 

-N [da\/jia)-piao)]pia) (1) 



where is the total number of particles in the system, P is the 
pressure, T is the temperature and the set {A/j} denotes the 
differences between fj{o), the chemical potential of a species 
with diameter a, and fj{Oo), the chemical potential of an (oth- 
erwise arbitrary) reference species: A/j(a) = fj{a) — /j(Oo). As 
we are dealing with hard-core particles, we choose our unit 
of energy to be equal to hgT. p{a) denotes the probability 
of finding a particle with diameter o. The set of thermody- 
namic fields {A/j} act as control parameters that determine 
the particle-size distribution. In the present work, we assume 
a quadratic dependence of A/j(o) on a — Gq: 



Ptu(o)-Moo)] = -(o-Oo)V2v 



(2) 



Where ^ = l/kBT. The parameter v determines the degree of 
polydispersity. At infinite dilution, the size distribution is di- 
rectly given by p{a) = cexp(— (a — Oo)^/2v). At finite con- 
centrations, the size distribution cannot be inferred directly 
from the functional form of A/j(a). Both the average parti- 
cle diameter and the actual polydispersity s (defined through 
= (o^) / (o)^ — 1) must be determined in the semigrand en- 
semble simulations. Once the functional form of A/j(o) has 
been specified, the semi-grand partition function S is a func- 
tion of N,P, r, V and Oq. 



l{N,P,T,v,ao)^ 
J dV J dr^ J da^ exp 



-^[PV + U{r^,a^)] 



(3) 
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The semigrand free energy Y is related to E through Y ~ 
—ksT InZ. To sample the configurations of the semi-grand 
ensemble, we use Metropolis-style Monte Carlo sampling of 
all variables that characterize a given configuration of the A^- 
particle system. In addition to the usual trial moves that at- 
tempt to change the particle coordinates {r^} and the system 
volume V, there are trial moves to change the diameter of a 
particle. As has been explained by Bolhuis and Kofke, it is 
computationally more efficient to combine volume-changing 
moves with particle resizing moves 1 3]. 

To calculate the chemical potential of the reference species, 
thermodynamic integration was used. As a reference state, we 
took the monodisperse hard-sphere crystal near coexistence. 
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for which the free energy per particle is accurately known fT] . 
In order to compute the change in free energy with P and v, 
we make use of the following thermodynamic relations: 



J N.T,<5Qy 



V 



Nldo'pia')^^^ (4) 



The semigrand free energy of an ideal, non-interacting system 
of polydisperse particles, is 



(O; - Oof 
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NkBT 

= N/Jidiao) = Gid -— ln(27i:v) 



(5) 



We can now employ the following scheme to compute /jgx (c^o) 
by thermodynamic integration, using as input our knowledge 
of the excess chemical potential /jgx.o of a monodisperse hard 
sphere system at pressure Pq: 



^ex(Oo) = A-ex.O + - dP' (v - ^ ^/ 



N Jq 



dV 



-Zi{Oi-Ooy- , NkBT\ 



2v'- 
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(6) 



realizable lattice in an orthorhombic simulation box; its value 
is well-defined by virtue of the extensivity of free energy. 

To calculate y^dd, we simulate a crystal with M lattice sites 
and M+ 1 particles, of which particle j has a scaled hard-core 
diameter aOj. The diameter scaling parameter a can be varied 
during the simulation, so that we sample the partition defined 
by 

f daZM.oAiM+\,PJ,ao,v,a) (8) 
Jo 

where Em,o,i(^+ i,PtT,V,Oo,a) is defined as in Eq.|3] but 
with configurational energy Uli^ ,0^^ ,aaj). We stress that 
particle j differs from the other particles only in the overlap 
criterion, not in the probability distribution that determines di- 
ameter sampling: for the overlap criterion, the particle radius 
of this particle is ac, whereas its weight in the Semigrand 
chemical potential distribution of Eq.|2]is still determined by 
o. 

During the simulation, we construct a histogram P{a\M + 
1,P,T,V): 



P{a\M+l,P,T,v) = 

Jq da'b{a - a')ZM.o.i{M + l,P,T,V,aQ,a) 

With this histogram we can calculate 



(9) 



B. Interstitial Concentration 

The methods that we used to calculate the concentration of 
point defects are similar to those discussed in Ref. 5. We first 
consider the free energy YM,nv,ni of a crystalline system con- 
taining M lattice sites, ny vacancies and rij interstitials. The 
total number of particles in this system is N = M + rij ~ ny ■ It 
is convenient to consider interstitials and vacancies separately. 

By analogy to the derivation of interstitial concentrations 
in monodisperse systems! 5], it is straightforward to show 
that the concentration of interstitials (xj) is given by xj w 
exp(— Py/), where yj is defined as yj — Ym,o,i — ^^m+i.o.o- It 
is convenient to rewrite yj as 



yi — iM,o,i 
= Ym. o.i 



■ iM+ 1,0,0 

■ iM.o.o + Ym, 0,0 — iM+ 1,0,0 



— iM,o,i — iM,o,o - 



A'id(c?o)+A'ex(<7o) 



= Ym.OA- iM,0,0 +A'id(Oo) -A'ex(C5o) 
= Jadd — jfexlC^o) 



(7) 



Here jadd is the free energy difference between a system 
with one interstitial and a perfect crystal plus one ideal (non- 
interacting) particle. The quantity yM+i,o,o, the free energy of 
a system with M + 1 lattice sites and no defects, is an abstract 
quantity that does not neccesarily correspond to a crystal with 



ygiow = -^B^ln 



f(a = l|M+l,P,r,v,Oo) 
P(fl = 0|M+l,P,r,v,Oo) 



(10) 



where jgrow is the reversible work needed to transform an in- 
teracting point particle (fl=0) into a particle with a hard-core 
diameter Oj (corresponding to «=!). In order to sample the 
full range of a-values from to 1, it is necessary to use bi- 
ased sampling. We employed multicanonical/umbrella sam- 
pling |6, 7] to generate P(fl|M,f,r,v,ao). 

To obtain the total interstitial free energy y^dd we must still 
add the free energy change associated with the transformation 
of a non-interacting particle into an interacting point particle. 
This free energy change is determined by the ratio of the vol- 
umes accessible to the two types of particles: 



Jadd Jgrow 



-kBT In 



(Vacc) 



V 



-kBTln{l-^) 



(11) 



where Vacc is the volume accessible to the point particle and 
r| denotes the volume fraction of the defect-free hard-sphere 
crystal. It is not necessary to confine the interstitial to a partic- 
ular Wigner-Seitz cell, as interstitials diffuse quickly through 
the system. If this were not the case, both the scaled and the 
unsealed particle would have to be confined to a particular 
Wigner-Seitz cell (or even, to one particular interstitial cav- 
ity). 
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C. Vacancy Concentration 

For the vacancies, we can get for the concentration xy ~ 
exp(-p3'v) (see Ref. 5), with yy = Ym+i^ifi - ^m.o.O and in- 
troduce the analogous free energy to jadd- 



yv 



iM+ 1,1,0 — Jm,o,o 

iM+1,1,0 — Jm+1,0,0 + iM+ 1,0,0 — iM,0,0 

Jm+1,1,0 - yM+1,0,0 + A'(<'o) 

iM+1,1,0 - Jm+1,0,0 +A'id(Oo) +A'ex(Oo) 



= - Y, 



M+ 1,0,0 - 



JM+l,l,0+A'id(<7o) ) +fe(Oo) 



-Jrem +A'ex(C5o) 



(12) 



In this case, jrem is the free energy difference between a 
perfect crystal and a crystal with one vacancy plus a non- 
interacting particle. 

If we assume that we can sample a system which can switch 
one particle between being a normal particle {b — b„) and a 
non-interacting particle {b = bi), we can introduce the equi- 
librium probabiUty P{b\M,P, r,v,Oo): 



-ksTln 
-ksT In 



P{b„\M,P,T,v,ao) 
Pibi\M,P,T,v,ao) 

{n{b„ ^ bi)) 



(13) 



where {K{bi b„)) is the mean transition probability from 
b = bi to b — b„. Because a real particle can always switch to 
a non-interacting, particle, we can reduce the expression for 

Jrem tO 



yiem = -kBT\\\{lz(bi -^b„)) 



(14) 



Now {K{bi b„)), the transition probability from a state of a 
system with a vacancy and a non-interacting particle to a per- 
fect crystal, is related to the probability Pin, for the insertion 
of a (normal poly disperse) particle into the vacancy: 



-kBTln{n{bi ^ b„)) = -kBT{lnP,„ 



(15) 



In practice, the simulation will consist of a collection of 
M — I normal particles and one ideal polydisperse particle 
which we keep in the Wigner-Seitz cell of the tracked vacancy. 
We then do multicanonical sampling, biasing on the number 
of overlaps that the ideal particle would create if it would be 
switched to a real particle, and get Pins from the probability 
to create zero overlaps. This scheme is essentially identical 
to that of Bennett and Alder j^, save for the multicanonical 
sampling. 
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FIG. 1 : Point defect concentration (x) versus polydispersity (s) 



5.8%. The latter value corresponds to the maximum polydis- 
persity attainable with the chemical potential difference func- 
tion used. Here, the polydispersity s is defined as the normal- 
ized second moment of the particle diameter distribution 



(CI) 



(16) 



All simulations were performed on 256(±1) particle sys- 
tems (a cubic fee 4x4x4 lattice); a simulation of a larger sys- 
tem in the monodisperse case in Ref.jSI shows that this particle 
number is sufficient for the required accuracy. For the (in- 
terstitial) calculation of ygmvi, the P{a\M +l,P,T) histograms 
were divided into 5 windows for which simulations were run 
in parallel. The multicanonical biasing weights were gener- 
ated starting with the weights for the monodisperse case and 
took 10-80 runs of4-105 MC sweeps (Monte Carlo cycles 
per particle) per CPU to converge. The final results were ob- 
tained using typically 80 runs of 4 • 10^ sweeps per CPU. In the 
case of vacancies there was one window for which about 20 
runs of 1 ■ 10*' sweeps were needed to equilibrate the weights 
after which about 40 runs of similar length were done for the 
final results. The equilibrium concentration of the two types 
of point vacancies as a function of different polydispersities is 
shown in Table U and Fig.fO 

The values of /Jex(oo), required for both the vacancy and in- 
terstitial concentration, were calculated using thermodynamic 
integration using the free energy differentials of Eq. |3 Inte- 
gration was done along the P-v points shown in Table |I] with 
20 steps between each step and 110^ averaging sweeps per 
step. 



III. RESULTS 

The simulations to calculate the point defect concentration 
were done at various points along the melting line of polydis- 
perse hard sphere crystals, as taken from Ref. 13, The points 
chosen give a polydispersity of approx. 1.5%, 3%, 5% and 



IV. DISCUSSION 

The simulation results show a dramatic increase in the in- 
terstitial concentration with increasing polydispersity, while 
the vacancy concentration remains roughly similar over the 
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0.54329 


0.54522(8) 


0.54641(6) 


0.55726(6) 


0.56997(6) 
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0.992 


0.967 


0.815 


0.589 


s 





0.015562(3) 


0.029974(7) 


0.05213(3) 


0.05755(5) 


A'ex 


17.071 


17.418 


18.308 


24.350 


37.516 


A'ex((cT)) 


17.1 
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17.8 
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22.5 
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7.92(1) 


8.098(9) 


8.77(2) 


13.68(4) 


26.1(2) 


^'grow 




32.2(1) 


30.8(2) 


29.5(2) 


40.5(1) 


xv 


1.10(2) •lO"'* 


9.55(9) -10-5 


8.3(2) 10-5 


4.6(2) -10-5 


5(1) -10-5 


XI 


2.7(4) ■ 10"*^ 


1.6(2) -lO^^ 


1.7(3) ■10"'' 


2.4(5)- 10"^ 


2.1(2)- 10^2 



TABLE I: Results for the vacancy and interstitial concentration for the polydisperse hard sphere system. The interstitial concentration for the 
monodisperse case was taken from Ref.?. All free energies are in units of and the pressure is in hgT /(5^, with the errors in the last digit(s) 
shown in brackets. Here, v is the polydispersity control parameter (see Eq.|2}, "n is the packing fraction, (o) is the mean packing fraction, i 
is the polydispersity, as defined in Eg. 1161 (o/) / (o) is the mean interstitial size relative to the mean particle size, /'ins is the particle insertion 
probability (see Ea. ll5> . jgrow is the free energy associated with growing an interstitial (see Eg. llOt . xy is the vacancy concentration and xj is 
the interstitial concentration. 
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FIG. 2: Normalized size distribution for total system and for the 
interstitials at polydispersities of 3.0% (left, v = 0.001) and 5.8% 
(right, V = 0.0056) 



full range of polydispersities. The increase in interstitial con- 
centration can be attributed to the size of the interstitials: if 
the particle size distribution has non-zero width, the intersti- 
tials are smaller than the mean particle size in the crystal, as 
is shown in Fig.|2] 

This size difference between interstitials and the surround- 
ing crystal is not an artifact of the simulation method: al- 
though the trial moves used in semigrand-canonical simula- 
tions are unphysical, the resulting size distribution of intersti- 
tials is real. The non-Gaussian particle size distribution in the 
crystal should be interpreted as a result of fractionation fsl 01 
of the coexisting fluid (with the same pressure and polydisper- 
sity control parameter v). 

The size distribution of the coexisting fluid is shown in 
Fig-E For small particle sizes, its value is slightly higher than 
the normal distribution, but at the peak of the interstitial size 
distribution, the difference in concentration is no more than 
7%. To a first approximation, the interstitial concentration in 
a crystal that has formed from a fluid in which the particle size 
distribution is exactly Gaussian, should be lower by the same 
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FIG. 3: Particle size distribution of the fluid coexisting with the high- 
est polydispersity solid (solid line, v = 0.0056, s — 0.058). The 
dashed line shows a normal distribution with the same first and sec- 
ond moment. The inset shows the probability distribution relative to 
the normal distribution. The vertical arrows mark the mean particle 
diameter of interstitials at the current polydispersity. 



amount. 

It must be noted, however, that once the crystalline phase 
starts occupying a sizable fraction of the system volume the 
size distribution will change and the interstitial concentration 
will probably be lower However, the exact size distribution in 
the crystalline phase is difficult to predict; the size distribution 
of the fluid itself will change as a result of the growth of the 
crystalline phase, and because of the high polydispersity of 
the coexisting fluid, the crystalline phase may be composed 
of several crystallites, each of which will have its own size 
distribution. 

The influence of the small particles on the interstitial con- 
centration can be illustrated by looking at the free energy of 
formation of a vacancy as a function of size. If we define a 
partial interstitial concentration x/(o), we can, as in Eq.0 ex- 
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FIG. 4: Interstitial free energy //(o) for different polydispersities i- 
as a function of renormalized particle size a/ (o). The v values for 
the different polydispersities can be found in Table Q The crosses 
denote the means of the interstitial sizes for the corresponding poly- 
dispersities. 

press it in terms of the free energy of formation //(o) and the 
chemical potential: 

x/(a) = exp(-p[//(o)-A(ex(a)]) (17) 

Assuming that the total interstitial concentration is the integral 
of the partial concentrations: 

XI = daxi{a) (18) 
Jo 

we can extract fi{o), the free energy associated with creating 
an interstitial of size a, because we know the chemical poten- 
tial distribution and the partial interstitial concentration. The 
values for //(a) at the polydispersities from TableUare shown 
in Fig 13 To be able to compare values of //(o) over a large 
range of a/ (a), the values for x/(a) in this figure were ob- 
tained by fitting the values from the simulations with locally 
skewed Gaussians 

x/(o) «flexp[-/7(o-(o/))2-c(o- (o/))3] (19) 

The fits work very well for the values of a which have been 
sampled during the simulation, and should yield meaningful 
results for the range shown in Fig. |4] 

The similarity in slopes and actual values of the //(o/ (o)) 
curves is striking; it means that, for the full range of poly- 
dispersities at which a crystal is stable, the partial interstitial 
concentration depends on the chemical potential distribution 
and an interstitial free energy which seems to be only weakly 
dependent on the polydispersity: 

with K = 74lkBT/al, ro = 0.338oo and/° = 1 1 .S/tgr as fitted 
parameters from the points in Fig0] Although the form of this 
equation was taken from the analytical estimate for the inter- 
stitial concentration of Ref. 5, which gives physical meanings 



to the values of K and ro and has reasonable agreement for ro, 
we stress that, here, K and ro are simply fit parameters. 

Because // (a/ (a)) hardly depends on the width and, pre- 
sumably, the shape of the particle size distribution, the small 
particle tail of the particle size distribution becomes crucial: 
those particles have the lowest // (a/ (a)) and will form the 
most important contribution to the interstitial concentration. 
For example, at the near-Gaussian polydispersity ofs — 5.2%, 
obtained by setting v ~ 0.004, practically all particles with di- 
ameter smaller than 75% of the mean particle radius are in- 
terstitials. This implies that the polydispersity, as measured 
by the second moment of the particle size distribution in the 
liquid, is not a good predictor for the interstitial concentration 
in the solid. The tail of the particle size distribution in the liq- 
uid is hard to measure, yet it is all-important for the interstitial 
concentration. 

In the case of vacancies, similar considerations apply in a 
slightly different form; the vacancy concentration depends on 
the chemical potential and the free energy of removing a par- 
ticle while keeping its lattice site. As argued above, they both 
stay relatively constant at melting for increasing polydisper- 
sities which causes the concentration of vacancies to remain 
roughly similar. 

To get an estimate for the interstitial concentration of a 
colloidal crystal in a suspension, the solely o-dependent ex- 
pression of Eq. |20| must be combined with an estimate for 
the chemical potential distribution /Jex(c?), which, in the more 
conventional ensembles of the experimental situation, does 
not only depend on the density and the mean particle size, 
but also on the subsequent moment of the particle size dis- 
tribution, the polydispersity |9, 10, 11]. An estimate for the 
absolute values of the chemical potential distribution can be 
obtained by combining Eq.|2]and the results of tableU 

In summary, we have shown that for polydisperse hard- 
sphere crystals along the melting curve, the interstitial con- 
centration increases dramatically (going up to 2%) while the 
vacancy concentration remains relatively constant. This can 
be attributed to the fact that, with increasing polydispersity, 
there is an increasing probability of finding a particle small 
enough to have an appreciable probability of fitting in a hole 
of the underlying crystalline lattice. 

This finding has practical implication for the preparation 
of colloidal crystals from slightly polydisperse solutions. As 
the presence of interstitials may affect the optical properties 
of colloidal crystals, it is important to control their concentra- 
tion. The present calculations show that the interstitial con- 
centration depends sensitively on the tail of the size distribu- 
tion in the liquid phase. Hence, the polydispersity as such 
does not provide a reliable criterion to predict interstitial con- 
centrations. Rather, it will be necessary to have an accurate 
representation of the functional form of the tail of the particle- 
size distribution (in particular, on the small-o side). 
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